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We present a computational study of the transport properties of campylotic (intrinsically curved) media. It 
is found that the relation between the flow through a campylotic media, consisting of randomly located 
curvature perturbations, and the average Ricci scalar of the system, exhibits two distinct functional 
expressions, depending on whether the typical spatial extent of the curvature perturbation lies above or 
below the critical value maximizing the overall scalar of curvature. Furthermore, the flow through such 
systems as a function of the number of curvature perturbations is found to present a sublinear behavior for 
large concentrations, due to the interference between curvature perturbations leading to an overall less 
curved space. We have also characterized the flux through such media as a function of the local Reynolds 
number and the scale of interaction between impurities. For the purpose of this study, we have also 
developed and validated a new lattice Boltzmann model. 

Many systems in Nature present either spatial curvature (e.g. curved space due to presence of stars 1 or flow 
on soap films 2 ) or geometric confinement constraining the degrees of freedom of particles moving on 
such media, e.g. solar photosphere 3 , flow between two rotating cylinders and spheres 4 " 6 , hemodynamics 
through deformable vessels 7 , fusion plasmas 8 , and flow through porous media 9 , to name but a few. In general, 
these systems force a fluid to move along non-straight trajectories (curved geodesies), leading to the upsurge of 
non-inertial forces. We will denote such systems as Campylotic, from the greek word koc/ukdXck; for curved, media. 
Due to the arbitrary trajectories that particles through a campylotic medium can take, depending on the com- 
plexity of the curved space, the flow through these media can present very unusual new transport properties. 
Campylotic media play a prominent role in all applications where metric curvature has a major impact on the flow 
structure and topology; biology, astrophysics and cosmology offering perhaps the most natural examples. Indeed, 
the flow through some specific campylotic media has already been studied, e.g. Taylor- Couette flow and flow 
through porous media. The Taylor-Couette flow was originally formulated for the case of two concentric, rotating 
cylinders 4,5 , and later extended to the case of spheres 610 . On the other hand, continuum descriptions of flow 
through homogeneous porous media have been already treated as changes in the fluid trajectories, by resorting to 
the concept of tortuosity, which is the ratio between the straight line distance to the curved path length (due to the 
presence of obstacles) between two points in the medium. The tortuosity is directly related to the permeability and 
porosity of the medium 9 . However, the flow through other complicated structures, like randomly located stars 
(which is an intrinsically curved space) or many biological systems, to the best of our knowledge, has never 
systematically been addressed before on quantitative grounds. 

Since, in general, this class of flows lacks analytical solutions, their study is inherently dependent on the 
availability of appropriate numerical methods. Flows in complex geometries, such as cars or airplanes, make a 
time-honored mainstream of computational fluid dynamics (CFD), a discipline which has made tremendous 
progress for the last decades 1112 . However, campylotic media set a major challenge even to the most sophisticated 
CFD methods, because the geometrical complexity is often such to command very high spatial accuracy to resolve 
the most acute metric and topological features of the flow. Therefore, in this work, we also present a new lattice 
kinetic scheme that can handle flows in virtually arbitrary complex manifolds in a very natural and elegant way, by 
resorting to a covariant formulation of the lattice Boltzmann (LB) kinetic equation in general coordinates and for 
curved spaces. 

As an additional feature, complex boundary conditions related with a specific geometry, e.g. surface of a sphere, 
or more sophisticated ones, like Mobius bands and the Klein bottle, can be treated exactly by cubic cells in the 
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contravariant coordinate frame, thereby avoiding the staircase 
approximations which would result from the use of cubic cells assoc- 
iated with euclidean geometry. For instance, in spherical coordinates, 
we can construct a lattice with cubic cells of volume dr X 36 X 3(j), 
being (r, 9, 0) the coordinates of each grid point. This feature can also 
be used to model the microscopic geometry of a porous medium by 
computing, analytically or numerically, the coordinate system that 
parametrizes the pore structure, and thus, to avoid the staircase 
approximation of the obstacles. The method is validated quantita- 
tively for very simple campylotic media with geometrical confine- 
ments by calculating the critical Reynolds number for the onset of the 
Taylor- Couette instability in concentric cylinders and spheres 5,610,13 , 
and applied to the case of two concentric tori. 

In this work, by using the new numerical scheme, we simulate the 
flow through campylotic media consisting of randomly distributed 
spatial curvature perturbations (see Fig. 1), in the laminar regime. 
The flow is characterized by the number of curvature perturbations 
and the average Ricci scalar of the space. The campylotic media 
explored in this work are static, in the sense that the metric tensor 
and curvature are prescribed at the outset once and for all, and do not 
evolve self- consistently with the flow. The latter case, which is a 
major mainstream of current numerical relativity 14 " 17 , especially in 
the turbulent regime, makes a very interesting subject for future 
extensions of this work. 

Results 

Kinetic theory in general manifolds. In order to study the campy- 
lotic media, we develop a lattice Boltzmann approach for curved 
spaces in general coordinates, taking into account the metric 
tensor and the Christ offel symbols V k -. The former characterizes 
the way to measure distances in space, while the latter is responsible 
for the non-inertial forces. The corresponding hydrodynamic 
equations can be obtained by replacing the partial derivatives by 
covariant ones, in both, the mass continuity and the momentum 
conservation equations. After some algebraic manipulations, the 
hydrodynamic equations read as follows: d t p + (pu% = 0, and 
d t (pu l ) + T l} - = 0, where the notation ;i denotes the covariant deriva- 
tive with respect to spatial component i (see SI). The energy tensor V j 
is given by, T ij = Pg^ + pu l ui - ^(^V / +/t^ / +^^ 5 where P is 

the hydrostatic pressure, u l the z-th contravariant component of 
the velocity, g the inverse of the metric tensor, p is the density of 
the fluid, and \i is the dynamic shear viscosity. 

Since lattice Boltzmann methods are based on kinetic theory, we 
construct our model by writing the Maxwell-Boltzmann distribution 
and the Boltzmann equation for general manifolds. The former takes 
the form 18 : 



/ eq =7^§7iexp 




Figure 1 | Streamlines of a three-dimensional fluid moving through a 
campylotic medium. The colors denote the scalar of curvature R' (blue and 
red for low and high values, respectively). The gray bubbles isosurfaces 
stand at 1/5 of the maximum curvature of the system. 
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where g is the determinant of the metric g t p and 6 is the normalized 
temperature. The macroscopic and microscopic velocities, u! and <f 
are both normalized with the speed of sound c s = \/k B To/m, k B being 
the Boltzmann constant, T 0 the typical temperature, and m the mass 
of the particles. Note that the metric tensor appears explicitly in the 
distribution function, due to the fact that the kinetic energy is a 
quadratic function of the velocity, u% = g^uK To recover the mac- 
roscopic fluid dynamic equations, we have to extract the moments 
from the equilibrium distribution function. The four first moments 
of the Maxwellian distribution function on a manifold are given by, 



p = jfd£ , pu^jfCdZ 



(2) 



(3) 



P 9(u^ k + iJg ik + u k g^ + pu i u^u k = jft'gfdl (4) 

These moments are sufficient to reproduce the mass and the 
momentum conservation equations. Here, for simplicity we have 
used d£ to denote d^d^ 2 d^ 3 and the Jacobian of the integration is 
already included in the Maxwell Boltzmann distribution, through the 
determinant term yfg. 

In the absence of external forces, in the standard theory of 
the Boltzmann equation, the single particle distribution function 
f(x\ £*, t) evolves, according to the equation, d t f + Qdjf = C(f), 
where C is the collision term, which, using the BGK approximation, 
can be written as, C= — (l/x)(f — f eq ), with the single relaxation 
time t. This equation can be obtained from a more general expres- 
sion, df /dt = C(f), where the total time derivative now includes 
a streaming term in velocity space due to external forces, 
df dx l dp 1 . 

-j- = d t f + -^-Cif+ —rVpif, wimp 1 the z-th contravariant compon- 
dt dt dt F 

ent of the momentum of the particles. Using the definition of velo- 
city, <f = cbf/dt, and due to the fact that the particles in our fluid move 
along geodesies, which implies the equation of motion 
dp l /dt= — T l kl p k p\ we can write the Boltzmann equation, using a 
procedure based on Ref. 19, as 



3/+a/-r;^V=e(/), 



(5) 



where we have used the definition of the momentum, p { = Note 
that the third term of the left hand side carries all the information on 
non-inertial forces. The Christoffel symbols and metric tensor are 
arbitrary and therefore we can model the fluid flow in curved spaces 
in general coordinates, whose metric tensor is very complicated and/ 
or only known numerically. 

Since the contravariant components of the velocity are free of 
space-dependent metric factors, they lend themselves to standard 
lattice Boltzmann discretization of velocity space. All the metric 
and non-inertial information is conveyed into the generalized local 
equilibria and forcing term, respectively. These features are key to the 
LB formulation in general manifolds. As an additional feature, even 
in flat spaces, complex boundary conditions related to a specific 
geometry, e.g. surface of sphere, in many cases, can be treated exactly 
by cubic cells in the contravariant coordinate frame, thereby avoiding 
stair- case approximations typical of cartesian grids. 

Studying the campylotic media. To study a genuinely campylotic 
medium, consisting of randomly located perturbations of spatial 
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intrinsic curvature, we define a coordinate system (x, y, z), such that 
its metric tensor deviates from the one of a flat space (Sy) in the form: 

gij = dij(^l — ao j exp ( — r n /r 0 ) ^ , where n labels each local 

curvature perturbation located at r n = (x n ,y n9 z n ), N is the total 
number of perturbations, r n = \r n \, and r 0 characterizes the size of 
the deformation. The choice of this metric tensor is to induce an 
intrinsic curvature in the space at each point f n . Note that the 
coefficient a 0 can be either signed, depending on whether a 
positive or negative curvature is imposed, respectively. In our 
study, we have chosen to work with positive values of a 0 , but its 
sign can be different according to the respective application. In 
order to understand the effects of the curvature, we will work in 
the laminar regime to avoid pressure losses due to turbulence. The 
effects of pressure drops in campylotic media due to fluid turbulence 
shall be addressed in future work. 

The Christoffel symbols are calculated numerically. The flux is 

calculated by the geometrical relation, <!> = J spu x \/g xx gdS, where 

S is the cross section at the location where the measurements are 
taken. Since the fluid dynamic equations only contain the metric 
tensor and its first derivatives (via the Christoffel symbol), it is 
natural to expect that the flow could be characterized by a quantity 
that contains the metric tensor and its first derivatives. Although 
the Christoffel symbols Y l - k meet this requirement, they are not com- 
ponents of a tensor, and therefore they are not invariant under a 
coordinate system transformation (physics should not depend 
on the choice of the coordinate system). An invariant, or tensor, 
that can be used to characterize the system is the Ricci tensor R ik . 
In this work, we use the Ricci scalar (curvature scalar) which can be 
calculated from the Ricci tensor, R t p by contraction of the indices, R' 
= g* j Rij. The relation between the metric tensor and Christoffel sym- 
bols and the Ricci tensor can be found in the SI. For sake of generality 
and to study our system using dimensionless units, we define a char- 
acteristic length X = ( V7iV) 1/3 , where V = L x L y L z is the total volume of 
our system, and the local Reynolds number Re = <S>r 0 /(tiS). 
Therefore, the dimensionless size of deformation is defined by 
e = r 0 /L We use a lattice size L x X L y X L z of 128 X 64 X 64, and 
t — 1. All quantities will be expressed in numerical units. To drive the 
fluid through the medium, we add an external force along the x- 
component, which in all simulations takes the value, f ext = 5 X 
10" 5 . The flux in flat space, i.e. in the absence of curvature perturba- 
tions is denoted by 

Shown in Fig. 1, are the velocity streamlines, the Ricci scalar R' and 
the high-curvature locations, represented by gray isosurfaces. Note 
that the streamlines are very complex, as the flow can orbit around 
the spheres before continuing its trajectory 20 ' 21 . This effect is due to 
the complex curvature of the space and has no relation to turbulent 
vortices, given the low Reynolds number (Re ~ 1). Also we can see 
how the curvature perturbations interact, creating non-spherical 
shaped isosurfaces. 

Fig. 2 shows the flux deficit <J> 0 - <I>, as function of the number of 
curvature perturbations, N We observe that the flux O decreases 
with N This effect is due to the interplay between the longer traject- 
ories that particles must take and the acceleration due to the non- 
inertial forces. Note that, in general, for systems with different 
configurations (e.g. negative a 0 ), we could expect that the combina- 
tion of the two effects might lead to higher flux by increasing N. We 
also see that the flux depends linearly on N for low concentration of 
curvature perturbations, and only sublinearly at higher concentra- 
tions. This is due to the fact that at low concentration, the average 
distance between curvature perturbations is large, and consequently 
each perturbation adds up as a single modification to the total spatial 
curvature. However, as the concentration is increased, the curvature 
perturbations start to interfere with each other and consequently the 
space becomes less curved (decrease of the overall Ricci curvature). 



e 
i 

c 

e 




Figure 2 | Flux through a campylotic medium consisting of randomly 
located curvature perturbations. Flux deficit O 0 — ® with respect to the 
flat case, as a function of the number of curvature perturbations for a 0 = 
0.01, r 0 = 2.0, and the varying between 0.1 and 0.4. The solid line is the 
analytical curve according to Eq. (6). Shown in the inset is the normalized 
average curvature scalar of the space, R/R max , and the normalized reduced 
flux 1 - Qq/® as a function of e (Re ~ 1 . . . 1000). Both Ricci scalar and flux 
deficit exhibit a maximum at intermediate values of e. Since the two 
maxima are slightly shifted with respect to each other, the reduced flow as a 
function of R exhibits two distinct functional expressions (see next Fig. 3). 



The flux is found to obey the following law, 

N/N 0 



(I) () (]) = A i 



1 + (N/N 0 ) 



(6) 



where A x = 163 ± 2 and N 0 = (1.54 ± 0.03) X 10 4 are fitting 
parameters. The parameter N 0 denotes a characteristic number of 
curvature perturbations, above which the sublinear behavior sets in 
(N > No). In dimensionless quantities, this is, 



(D () (D = Ai 



l + (6/e 0 ) 6 



where e 0 = 0.616 + 0.005 is the characteristic value above which the 
sublinear behavior starts to be dominant, and provides an "effective 
radius", r 0 /£~3.25, for each curvature perturbation. 

In the inset of Fig. 2, we observe that by fixing the number of 
curvature perturbations N = 1024 and the strength a 0 = 2 X 10" 5 , 
and changing the range of the perturbation, e with 2 = 8, the differ- 
ence ®o - O presents a maximum for a given e = e c ~ 4. Note that this 
value is similar to e 0 . Due to the fact that by increasing r 0 the flux O 
decreases, one could think that the Reynolds number would achieve a 
plateau. However, by analysing our numerical data we have not 
found such behavior. Furthermore, another interesting result is that 
the average dimensionless curvature, here defined as R = — 10 6 1 2 < 
R' > (where < ... > means average over space), shows the same 
qualitative behavior. Since by increasing e the metric tensor compo- 
nents decrease monotonically, this maximum is due to the 
Christoffel symbols (or non-inertial forces), which can be character- 
ized via R. However, the maxima are slightly shifted, due to the fact 
that the Ricci scalar does not uniquely determine the metric tensor 
and Christoffel symbols, the quantities that play a key role in the fluid 
dynamic equations. Taking into account this effect, we can plot the 
flux deficit <E 0 - ® as a function of R, and find that, indeed, for e < e c , 
the flux decreases by increasing the average curvature R with a dif- 
ferent law than for the case of large values of e > e c (see inset of Fig. 3). 
This gives rise to a loop-shaped curve, the reason for this behavior 
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being that the metric tensor is different for e < e c and e > e c , even if R 
takes the same value. However, in both cases, the system shows that 
higher values of the average curvature R always result in a lower flux. 
The behavior of the flux for e < e c is well represented by the following 
law: 



(D () (I) = A 2 - ( 1+ — ), 

\ Roy 



and for the case of e > e c , 



(]) () cD = A : 



R 



Ro+R 



(8) 



(9) 



where R 0 = 5.2 ± 0.1, A 2 = 50 ± 2, A 3 = 154 ± 4, and ®' = 5 ± 1. 
The quantity £0 is related to the maximum curvature achieved by the 
system and the intersection of the two laws (see Fig. 3). The other 
interesting quantity is O', which represents the difference of flux 
between £^>e c and e^Ce c , when the curvature scalar becomes zero, 
and it is due to the fact that in both cases, although the space has no 
curvature, it has nonetheless different metric tensors. 

Finally, in order to give a more general description to our campy- 
lotic medium, we perform a detailed study of the the flux <I>, as a 
function of two dimensionless number, Re and e. Since Re is not the 
only dimensionless number that characterizes the flow, we can 
expect different behaviors depending on how the Re is varied. 
Therefore, we will change Re by modifying both, the perturbation 
range r 0 and the viscosity fi, separately. In Fig. 4, we can appreciate 
the dependence of the flux deficit, <X> 0 — O, as a function of e for a fix 
Reynolds number, Re = 3. Note that fixing the Reynolds number and 
changing the ratio e, implies to vary the number of impurities per- 
turbations. As a consequence, in order to achieve a broad range of 
values, we decrease the amplitude a 0 , taking the value of 2 X 10~ 6 . In 
the domain ee[l, 10], we found that the flux decreases following a 
power law with exponent 3. 

In the inset of Fig. 4, we report the behavior of the flux deficit as a 
function of the Reynolds number by changing r 0 , for e = 5 and dif- 
ferent values of the fluid viscosity fi. Here, in order to vary Re, via r 0 , 
for fixed e and fi, we need to modify both, N and r 0 . In the regime of 
Reynolds numbers studied in this work, and for a fixed viscosity, the 
flux increases with the Reynolds number almost linearly upto some 
critical number, Re c , above which a sublinear behavior is found 



nearing a plateau. The flux deficit obeys the law, 



l + (Re/Re c ) 



(10) 



with3 2 = 0.259 ± 0.006 and£e c = 5.8 ± 0.1 for fi = 0.18, B 2 = 0.289 
± 0.008 and Re c = 0.75 ± 0.02 for fi = 0.55, and B 2 = 0.16 ± 0.005 
and Re c = 35.8 ± 0.9 for fi = 0.07. Note that Eq. (10) is satisfied 
regardless the viscosity of the fluid, and we can conclude that by 
decreasing the viscosity, we increase the critical Reynolds number 
at which the sublinear regime begins. On the other hand, the insens- 
itivity of the flux deficit at large Reynolds numbers (by changing the 
interaction range, r 0 , and therefore, large values of r 0 ) is due to the 
fact that the average curvature of the system decreases, leading to a 
flat space, and O (1 + 5 2 )<I>o- 

Discussion 

Summarizing, we have explored the laws that rule the flow through 
campylotic media consisting of randomly distributed curvature per- 
turbations, and shown that, for the configurations studied in this 
work, curved spaces invariably support less flux than flat spaces. 
Furthermore, the flux can be characterized by the Ricci scalar, a 
geometrical invariant that contains the metric tensor and 
Christoffel symbols. The trajectories of the flow can become very 
complicated, due to the total curvature of the medium, presenting, 
in some cases, orbits winding several times around regions with high 
curvature. To add generality to our study, we have also analyzed the 
flux transport across the campylotic medium as a function of dimen- 
sionless numbers, e and Re. Furthermore, we have found that the 
different laws that characterize the campylotic medium are valid 
regardless of the viscosity of the fluid, as far as the laminar regime 
holds. Indeed, further extensions to consider turbulent flows will be a 
subject of future research. 

To calculate the flux in campylotic media, we have developed a 
new lattice Boltzmann model to simulate fluid dynamics in curved 
manifolds using an arbitrary coordinate system. The model has been 
successfully validated (see SI) on the Taylor- Couette instability for 
the case of two concentric cylinders and spheres, the inner rotating 
with a given speed and the outer being fixed. We also studied the 
Taylor- Couette instability in two concentric rotating tori, finding 




Figure 3 | Flux deficit, <£ 0 - ^> as a function of the average curvature, R, 
for large and small values of e. We have fixed Oq = 0.00002, N = 1024 (Re 
~ 1 . . . 1000). The solid lines denote the analytical curves according to Eqs. 
(8) and (9). The inset shows the loop which arises by parametrizing the 
flux- curvature relation in terms of e. Here, e c is the radius at which the Ricci 
curvature attains its maximum upon increasing e. The lower and upper 
branches correspond to e<e c and e>e c , respectively. 




Figure 4 | Flux as a function of the dimensionless numbers Re and e. Flux 
deficit, <X> 0 — O, as a function of the dimensionless number e for Re = 3, 
which decreases following a power law with exponent 3, in the studied 
regime. In the inset, we can observe the dependence on Re by fixing e = 5 
and varying r 0 , for different viscosities, v = 0.55 (green squares), 0.18 (blue 
circles), and 0.07 (red diamonds). 
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that the critical Reynolds number for the onset of the instability is 
about ten percent larger than the one for the cylinder. By solving the 
Navier-Stokes equations for curved spaces in contravariant coordi- 
nates, which can be represented on a cubic lattice precisely in the 
format requested by the lattice Boltzmann formulation, the present 
model opens up the possibility to study fluid dynamics in complex 
manifolds by retaining the outstanding simplicity and computational 
efficiency of the standard lattice Boltzmann method in cartesian 
coordinates. 

Methods 

In order to formulate a corresponding lattice Boltzmann model, we implement an 
expansion of the Maxwell-Boltzmann distribution, Eq. (1), in Hermite polynomials, 
so as to recover the moments of the distribution function up to third order in 
velocities, as it is needed to correctly reproduce the dissipation term in the hydro- 
dynamic equations. The expansion of the Maxwell-Boltzmann distribution was 
introduced by Grad in his 13 moment system 22 . Since this expansion is performed in 
velocity space, and the metric only depends on the spatial coordinates, we expect such 
an expansion to preserve its validity also in the case of a general manifold. We have 
followed a similar procedure as the one described in Refs. 23, 24. 

For the discretization of the Maxwell Boltzmann distribution (1) and the 
Boltzmann equation (5), we need a discrete velocity configuration supporting the 
expansion up to third order in Hermite polynomials. Our scheme is based on the 
D3Q41 (three dimensions and 41 velocity vectors) lattice proposed in Ref. 25, which 
corresponds to the minimum configuration supporting third-order isotropy in three 
spatial dimensions, along with a H-theorem for future entropic extensions 26 of the 
present work. 

In the following, we shall use the notation c\ to denote the z'-th contravariant 
component of the vector numbered X. Thus, the discrete Boltzmann equation for our 

model takes the form,/;. (V + c\dt,t+ 5t) -fx (x\t) = -^(f k -/ ; eq ) + 5tTx> where 

T i is the forcing term, which contains the Christoffel symbols, and/ eq is the discrete 
form of the Maxwell-Boltzmann distribution, Eq. (1). The relevant physical 
information about the fluid and the geometry of the system is contained in these two 
terms. The macroscopic variables are obtained according to the relations, 

^= EL/^ h - EL-^- 

The cell configuration D3Q41 has the discrete velocity vectors: (0, 0, 0), (±1,0, 0), 
(±1, ±1,0), (±1, ±1, ±1),(±3,0, 0), (0, ±3,0),(0, 0, ±3), and (±3, ±3, ±3). The 
speed of sound for this configuration is c] = 1 — a/2/5. With this setup, and taking 
into account that the vectors <f and u l are normalized by the speed of sound, we obtain 
the following discrete equilibrium distribution, 



fT 



5 c\ u 1 
2 +2 t 



uu 
' 2cf 



i44 , \ (<W) 

2 c? 2 cf 



2 2c] ' 6c 6 s 

(C;W ! )(W ; V) (4 M ') 



2 cf 



2ct 



(4/4- 



(11) 



where the weights w k are defined as, W( 0 ,o,o) = [5045 — 1507^/10 J , 

37 91 If r-\ 1 / r- \ 

Whoo) = — 7= > 10) = — 55-17V10 J, W(i 1 d = 233V10-730 , 

{ " } 5y/l0 40 ( " j 50 V / {W) 1600 V / 

w ( 3 >°>°) = T^n" (295-92^), and w {w) = (l30-41^To) . The forcing 



16200 
term takes the form 



2c? 



c\F[ \ 



2c] 
(c\u l )c\F\ 



2c] 



u l F\ 



(12) 



with F\ = — Tj k ^^l and 3 kl the Kronecker delta. In the presence of an external force 
F ext , this simply extends to F\-^F\ + F l exV 

In order to recover the correct macroscopic fluid equations, via a Chapman- 
Enskog expansion, the other moments, Eqs. (2), (3), and (4), also need to be repro- 
duced. A straightforward calculation shows that the equilibrium distribution function 
fP meets the requirement. The shear viscosity of the fluid can also be calculated as 
fi = p(x — l/2)c 2 s St. In this way one can calculate the fluid motion in spaces with 
arbitrary local curvatures. 

We have measured the efficiency vs. standard LB and the resulting overhead (about 
3) is almost entirely to be ascribed to the fact that, by relativistic necessity (third order 
moments to be matched), we work with 41 velocities. Although a detailed head-on 
comparison remains to be done in future work, such overhead appears perfectly 



acceptable, especially in view of future applications to relativistic dissipative hydro- 
dynamics in highly complex manifolds. More details about the discretization of this 
model on a lattice and the numerical validation can be found in the SI. 

All simulations in this paper are performed using the smallest system size that 
keeps the physical results unchanged at increasing the grid resolution. The reason for 
this choice is that we need to calculate the metric tensor and the Christoffel symbols 
locally, which calls for a computational compromise. On the one hand, storing the 
metric tensor and Christoffel symbols as arrays, minimizes the CPU time, at the price 
of increasing memory requirements. On the other hand, by computing these 
quantities "on the fly", memory request is significantly reduced, to the detriment of 
computational time. For the simulations presented in this work, we have used the 
former approach (metric and Christoffel symbols as stored as arrays), so that we 
cannot compute very large system sizes in a reasonable computational time. However, 
finding the optimal tradeoff between the two approaches above will be a subject of 
future work. 
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